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A  Report  on  the  rJuinerlcal  Solution 
of  the  Geostrophic  Conservation  Equation 

David  A.  Levins 

I.   Introduction 

At  present  i/jork  is  progressing  under  the  guidance  of 
E.  Isaacson  and  G.  K.  horikawa,  and  others,  aiming  at  a  feasible 
numerical  method  of  solvin,c-  the  geostrophic  conservation  equa- 
tion.  It  is  believed  that  solutions  of  this  equation  may  con- 
tribute to  the  understanding  of  the  problem  of  meteorological 
forecasting,  especially  when  applied  to  the  problem  of  predicting 
the  path  of  large  disturbances  such  as  hurricanes. 

In  this  paper  we  shall  deal  mainly  with  describing  the 
numerical  methods  which  \-ient    into  solving  the  equations.   Other 

papers  presently  in  preparation  will  discuss  the  meteorological 

9 
implications  of  results  obtained  so  far. 

In  the  folloX'.jing  we  shall  indicate  the  physical  model  used 
and  give  a  statement  of  the  derivation  of  the  basic  equations. 
For  a  detailed  treatment  of  this  derivation  one  should  consult 
the  paper  of  reference. 

The  physical  model  used  is  that  of  a  single  layer  of  in- 
compressible, inviscid,  homo;~eneous  atmosphere  at  constant  den- 
sity, with  a  free  surface  at  a  height  which  includes  most  of  the 
atmosphere  when  equivalently  reduced  to  constant  density.   The 
height  of  the  free  surface  is  small  compared  l^Jith  the  average 
wave  length  of  large  scale  disturbances,  and  hence  the  shallow 


■  i  :f- 
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water,  gravity  wava  approximations  may  be  used.     Alternatively 
one  may  use,  if  necessary,  a  model  of  several  layers  of  constant 
densities  with  a  free  surface  at  the  topmost  layer  and  inter- 
faces between  layers.   In  particular  we  also  use  the  two  layer 
model  with  a  deep  layer  of  less  dense  warm  fluid  overlying  a 
relatively  thin  layer  of  denser  cold  fluid. 

Since  the  phenomena  under  study  are  local  compared  with 
the  surface  area  of  the  earth,  we  take  instead  of  a  spherical 
ground  surface  a  plane  tanr/ent  to  the  earth,  with  the  point  of 
tangency  centered  at  the  locale  of  interest. 

The  derivation  of  the  basic  equations  is  from  the  shallow 
water  gravity  wave  equations.    An  exact  solution  of  these 
equations  is  given  by  a  free  surface  of  constant  height  and  an 
atmosphere  at  rest.   Because  of  the  slow  movement  of  large  scale 
atmospheric  phenomena  this  basic  solution  is  appropriate  for  a 
perturbation  expansion.   The  quantities  used  here  are  the  first 
order  terms  of  a  perturbation  on  this  state  of  rest.   The  order 
of  the  terms  in  this  expansion  is  selected  on  a  strictly  mathe- 
matical basis  which  is  implied  by  the  conservation  of  certain 
quantities  and  the  assuinpt  io]i  of  peostrophic  flow.   The  result- 
ing equations  for  the  first  order  terms  of  the  perturbation  ex- 
pansion are  a  nonlinear  set  and  are  called,  collectively,  the 

seostrophic  conservation  equation.   These  equations  are  related 
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to  those  recently  considered  by  B.  Bolin.  " 


II.   statement  of  the  Mathematical  Problem 

Set  up  a  rectangul'ar  coordinate  system  in  the  plane  of 
tangency,  with  the  oririn  at  the  point  of  tanp-ency.   Let  the  x 
coordinate  increase  eastward  and  the  y  coordinate  increase  north- 
ward.  Let  the  time  be  t.   Set  T  =  at,  where  g  is  a  scaling 
constant.   Let  h^  '{x,j,X)    be  the  height  to  first  order  of  the 
free  surface  above  or  below  a  constant  plane  surface  of  height 
h^.   Here  h^  is  the  constant  height  of  the  free  surface  for  the 
flow  at  rest.   The  quantity  \\i^    ■    is  defined  by 

where  g  Is  the  acceleration  due  to  gravity  and  f  is  the  Corlolls 
parameter: 

(2)  f  =  2Alsin  f   . 

Here  -H.  is  the  angular  velocity  of  the  earth  and  ^  is  the  lati- 
tude.  For  convenience,  we  choose  f  as  sensibly  constant,  eva- 
luating it  at  some  central  point  of  interest  —  this  is  not 
essential  to  the  aiethod  of  solution,  as  it  may  easily  be  allowed 
to  vary. 

We  now  wish  to  study  the  atmospheric  flow  for  times  T, 
0  <  X  <  T  at  a  level  high  enough  to  exclude  the  turbulence 
effects  of  the  ground  —  say  a  height  corresponding  to  500  milli- 
bars of  pressure  (approximately  19  to  20  thousand  feet).   Let 
the  com.ponents  of  the  first-order  horizontal  velocity  at  this 
level  be  u   '^(x,y,t)  and  v    (x,y,t).   The  basic  equations  of 
flow  are : 


(3)  u^l)  =  -t^^ 

-1 

Here  letter  subscripts  denote  partial  differontiation,  and  we 
have : 

a)  V^  1^^^   =   ?^^^  +  ?^^^ 
—XX     --yy 

b,  S^=(  L.u<^)(  ),.v(^)(  ,^ 

2  2    f 

c)   Vf   is  a  constant.   V^;  =  -t~- 

The  letters  f,  g,  and  h  e.re    as  before.   VJe  now  drop  the  super- 
script (1)  for  convenience  and  speak  only  of  the  first  order 
terms  of  the  perturbation  expansion.   The  systom  (3)  is  to  be 
solved  under  the  followin;'"  conditions: 

(ij.)   Domain:  The  domain,  D,  of  the  solution  is  the 

rectanple  x  J^x<x.,  ;y  <y<yj  for 
times  0  <  X  <  T. 

Initial  Values:    I'^jy.x)  is  given  as  a  function  of  (x,y) 

in  D  for  T  =  0. 

Boundary   Values:      j_{x,y,T.)    is    given   on   the   boundary    for   all 

T.  in   0   <  t    <   T.      \/  J{x,j ,i)    is    given   on 
the    boundary    for  all   T,    only    at    those 
points    of   the   boixndary    for  which,    at 


time  t,  the  velocity  is  directed  Into  the 
domain, D.   (Alternately,  one  could  pre- 
scribe  "^  1  -  V  ?  in  the  same  manner, 
since  ^   is  pi  escribed  on  the  boundary.) 

With  these  mixed  initial  and  boundary  values  prescribed 

the  solution  of  system  (3)  has  been  shovjn  to  exist  and  is 
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unique.     The  restriction  specifying  the  vorticity,  "xZ  1,  on 

the  boundary  only  at  times  and  points  where  the  fluid  enters 
the  domain  comes  from  the  following  consideration.   Since  _f  is 
known  at  all  tim.es  on  the  boundary,  the  quantity  "V  f  -  <  f  is 
known  on  the  boundary  whenever  V  V  is  known  there.   Hence  spe- 
cifying  \7  ?  is  equivalent  to  specifying  "v^  ?  -  h  !•   But  equa- 
tions  (3)  state  that  \7  _|  -  ^  ^  is  at  time  X,  constant  along  any 
particle  trajectory.   At  time  x    if  \7  Y  -  k  jj  is  specified  on 
the  boundary  at  a  point  P  where  the  flovi  is  out?,oin.y,-  this  spe- 
cified value  K'ould,  in  general,  contradict  that  computed  in  the 

interior  of  D  glon.'r  the  particle  trajectory  leaving  the  domain 

p     '^ 
at  P.   On  the  other  hani,  at  time  t,  \7   }[  -  k.""?  mu3t  be  speci- 
fied at  a  point  F  oxi   the  boundary  wheie  the  flow  is  incoming, 
since  it  cannot  be  computed  from  the  interior  —  in  fact,  the 
specified  value  of  \/'~^   -  K.  J[  at  P  determines  the  value  of 
\7  _L  "  ^'-'"^  ^'t  every  interior  point  included  in  the  particle 
trajectorj?-  through  P. 

For  the  purpose  of  studying  large  scale  centered  disturb- 
ances such  as  hurricanes  we  divide  the  quantity  ^   into  two  parts: 

(5)  I(x,y,r)  =  I^(x,y,T)  +  l^(x,y,T)   . 


Here  ^^  is  an  unknown  function  to  be  solved  for  and  represents 
the  stream  function  of  the  regular  or  steering  flow.   ?  is 
derived  from  the  atmospheric  vortex  theory  as  the  known  function: 

(6)  I,(x,y,t)  =gK^(Kr^)   . 

K   is  the  second  solution  of  the  modified  Bessel's  equation.   It 

is  a  positive  function  with  a  logarithmic  singularity  at  the 

origin  and  dies  out  exponentially  for  large  argxmient.   y  is  si 

positive  constant,  the  strength  of  the  vortex,  given  in  reference 

8.   And  r   is  defined  as: 
o 


(7)    a.        r^  %/^Xo)'  ^  (^-yo)'   * 

Here  x  =  ^'q^"^)  ^^^  y  ~  y  ^'^"^    ^^'^   ^^^®  coordinates  of  the  center 

of  the  vortex,  which,  from  the  atmospheric  vortex  theory,  moves 
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according  to  the  ordinary  differential  equations: 

dx 


(7)   b. 


^  -~   -i\}^ 


dy 

— --  =  ( ^f  ) 
d'i.     ^-^I'x 


^o'^C'^ 


o'  o'  ^ 


It  is  clear  from  the  above  that  the  motion  of  the  vortex 
is  determined  by  the  stream  function  of  the  regular  or  steering 
floiij,  jn  •   However,  it  must  be  pointed  out  that  the  vortex  also 
modifies  the  steering  flow.   The  function  ^     is  separated  out 
from  the  overall  flow  _t  ^s  a  known  function  for  the  purpose  of 
keeping  track  of  it  in  calculations  which  would  otherwise  wash 
it  out.   This  will  become  clear  in  the  sequel  when  it  is  realized 
that  in  practical  computation  the  spacing  between  mesh  points 


represent  distances  on  the  order  of  200  miles  —  the  center  of 
the  vortex  must  be  located  much  mora  accurately  than  this  mesh 
si7.e  indicates. 


III.   The  Numerical  Method  of  Solution 

For  the  purposes  of  computation  we  write  equation  (3)  as 
follows 

TJ^   .  y,2|  .  P(.x,y,t) 
dP 


(8) 


d'C 


=  0 


riow   since   ?  =   f^   +  I.    and  Jq  =  ""^  '^o^'^^o^    ^^    ^   solution   of 
V^I  -  k^i  =    5(x,y),    we   have  V%   -   V?l^  =   0  for  x  ^  x^,   y  /  y^, 
Hence   we  may  write    sj'-stem   (3)    as: 


(9) 


2tc 


..    V%  -  y'.%  =  p 


b.      P      -?P     +?P     =0 
t        -^y  X       ^x  y 


c.  y 

— o 


dx 


J-   K    {\\.T    ) 
2'ii      o        o ' 


dy 


d-t 


-di). 


(Il)x 


X     ,  V     , ' 

o'-  o' 


X     ,  Y     ,  T 
O  '  "  o ' 


Note   that    in   equation    (9b)    the   velocities    u  =    -f    , 


— X 


are  the  velocities  of  the  total  flow,  i.e.  f^  =    [(I^)^  +  ^^o^x^ 
and  ^  =  [(I-i)   +  (I  ),J-   If  it  were  not  for  these  Quantities 
deriving  from  the  total  stream  function  we  could  solve  the  equa- 
tions first  for  the  regular  or  steerin?  flow  and  determine  the 
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I 
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trajectory  of  the  vortex  as  an  afterthought.   But  this  is  not 

possible,  since  in  (9b)  the  velocity  field  due  to  the  vortex 

enters  also  iuto  the  computation  of  the  steering  flow. 

V'e  now  set  up  a  grid  mesh  in  the  rectangular  domain  D  and 

replace  the  continuous  variables  (x,y)  by  (x.,y.)  as  follows: 

J  "^ 

(10)   a.   X.  -  x^  =  j/^x      ;        j  =  0,1,2, ... ,M+1 
J     o 


X 


-  ^.  =  (^'+i)A 


w 


d   "c 

b.   y^  -  Yi  =  iAy   ;        i  =  0,1,2,...,N+1 

^d  -  ^c  =  (^'^'^i^Ay 

c-       'tjc  "  ^^^'^      '        ^^   "  0,1,2,...,T 

ith  /\x,    A. y ,  /\'l,  fixed  consto.nts.   "e  now  write  the  functions 


of  continuous  variables  as  functions  at  the  discrete  points  of 
the  mesh,  i.e.,  we  write  f(x.,y^,ljj  as  f  .  ^.   The  only  quantity 
we  leave  as  a  function  of  continuoiis  variables  is  the  coordinates 
of  the  position  of  the  vortex  at  discrete  time  steps  k  —  which 

we  write  as  (^qsYq)  • 

The  method  of  continuing  information  from  one  time  step 
to  the  next  is  the  following.   The  function  ij^)  .   .    is  determined 
only  from  (x  ,y  )  .   Initially  we  have  given  ^-  4  and  the  posi- 
tion of  the  vortex  i^^}!;^)    •   Hence  we  can  compute 
(?  )°  .  =  IjJ^  .  -  (^  )9  ..   Then  we  can  compute  f9  ,  from 
^2^m  ^0  _  _  ^2^w  )0  ^  ^  pO  _^   2^  ^^^^^  ^^^  ^^^   assume  we  have 

initially  eiven  P.  .,  (1^  )  ,•  ,•  and  (x  ,y  )  .   Wow  suppose  we  have 

given  at  time  step  k,  the  quantities  F.  .,  (^-,  )  •  .,    (^  ,y  )'. 

Jj-'-     -'-Jj-'-     '-'   '-' 


1.  We  first  compute  (f„).  .  from  knowledc^e  of  (x  ,y  )  ,  i.e., 

r     ^  -1-0  '  j ,  1  '      ^  o '  '^  o   '      ' 

^o  ^  "si  ^"^o^^'^o^*   ^~'®®  section  VII). 

2.  Ve   then  compute  ?^  .  =  (][,  )^  .  +  (f  )^  .  and  its  derivatives 
(j.-,J,-  ^  >  (2,.)  .,•  ,•  usinp:  central  differences. 

3.  Now  using  equation  (9b)  and  the  known  boundary  values  of 
P'^  .  for  incoming;  flew,  we  find  F.  ..   (See  section  IV). 

L|..   Using  F.  .,  the  known  boimdary  values  of  (j)  .  .  and 
J  >  J-  J  J 1 

equation  (9a)  we  solve  for  (Y-,)"'^  ..   (3ee  section  V). 

5.   We  now  find  (x^.y^)^'^  by  using  (l^^^i,  ^i^^i^i  "^^^  ^^^o'^o^' 
( See  section  VI ) . 

Hence  after  step  5,  we  have  F'.  .,  (?-|  )  .  •  and  (x  ,y  ) 
^''e  can  then  continue  to  the  next  time  step. 

I V .   Soluti on  of  the  Hyperbolic  F quation  (^b) 

(9b)  F^  +  uF   +  vP  =  0   . 

•      X     y 

W©   replace   the    time   derivative   b^"-   a   forward   difference    and 

the    space  deriv^.tives   by  a    trackward  difference    or  a   forward  dif- 

■i 

ference  depending  on  the  s^'gn  of  the  velocity  com.ponent. 

That  is: 


(11)   a.   (FV)^  .  =  -^  (F-;^^  -  ?^^  .) 


b.   u^  .(F  )^.    .    '-   -^    lu'^  .|(^^  ,  -  F^        ,      ) 
J'^   '^  J'^   Ax   ^'^        -i    (j-sgn  u\^  .),i 


(11)       c.      v^    .(F    )^    .    =   J-    |v^    .|(P'^    ,    -    F^ 

J'^     y    J'^        /xy        ^'^        ^'^  j,(i+3.n  v^^    .) 

J  >  ^ 


10 
) 


d.      sgn  f  =    < 


1    for  f   >  0 


V 


1   for   f  <  0      . 


v+i 
Hence    solving  for   P".    .    we  have 


(12)      P 


k+1 


1     .    _ |u 


Ay 


V   .      , 


-,K 


j,i 


^  A^  ,  ic     I  ^k  _^  At  ,  k     I  ^k 

A^        J,i        Vj-sgnu),i        ^,y        j,i'       j,(i+sgn  v) 

Now  assume,  only  for  purposes  of  estiraating  the  error, 

k        k 
that  the  coefficients  u.  .  and  v".  .  are  known  functions  and 

hence  the  equations  are  linear.   We  let  (P.  .  -   F.    .)    =    e.    ., 

J,i     J,i      J,i 

where  P.  .  is  the  true  solution  of  the  differential  equation. 
Then 

(1.3) 


,k+l 


Ai    , 
1 |u 


Ax 


A 


J.- 


Zl 


A^ 


j-sf-n  u)  ,i 


At 


V      e 


j, (i+sgn  v) 


o(At  )  +  o(A"^A^)  -«■  ^'(A^Ay)    • 


Nox^!  we  note  that  the  error  at  the  (k+l)st  step  is  the 

weighted  averarre  of  the  error  at  three  points  of  the  k-th  step 

with  weirhts  that  sum  to  unity.  If  each  of  the  tjeights  is  non- 
ner:atlve  we  can  conclude: 

(11+)   max  le^""^!  <  max  |  e^  .|  <•  O(Al^)  +  0{/^x/\x)    +  O(AxAy) 


J,l 


J5I 
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Now  let 


h:       I  k   I 
e  =  ma  X   e  .  . 

Then  itorating  equation  (lii)  k  times  ve  get: 

(15)   e'^  <  e°  +  t[0(/\0  »-  0{/^x)    +  0(/\y)],   since  k/\T  =  T.  . 

This  gives  an  estimate  for  the  erroi'  at  tirae  step  k  in 
terms  of  the  initial  round-off  error  and  the  truncation  error. 
We  note  th^t  if  the  initial  round-off  error  vanishes  then  the 
error  at  time  r  tends  to  zero  with  {/\x,    /\y ,  /\t,)  . 

The  condition  that  the  weights  in  equation  (13)  be  posi- 
tive is  satisfied  bv  the  computational  stability  condition: 


r 


/    Ax       /V  /7"k  .Z   ^    /    k      sZ 

V  ^-^  ^— ^       max   /  ^  u .  .  )   +  ( V  .  .  ) 


(16)    /  i  —  f   +  (r___) 


j,k^      J.l  J,l 


^ •   The  Solution  of  the  Elliptic  Equation  (9a) 

(9a)  V%  -  n\   =   F      . 

Assume  we  are  at  time  X.   For  convenience  we  now  reserve 
the  superscript  for  other  purposes  and  drop  the  subscript,  1, 
from  ^-.  .      For  the  solution  of  this  equation  we  take  the  finite 
difference  scheme: 


12 

whe  re 

^■'hen  ^.    .  is  evaluated  on  the  first  or  last  row  or 
J  >  - 

column  we  take  the  known  boundary  values  neighboring  j .  .  to 

J  >  •'- 

the  right  band  side  of  equation  (17)  and  denote  the  right  hand 
side  as  G.  ..   Then  we  may  write  this  system  of  simultaneous 
linear  equations  in  the  following  manner.   Let  each  colunin  of 
ISi  values  of  ? .  .  be  denoted  by 


and  the  column  of  K  such  successive  columns  by 


/ 


Similarly  we  set 


and 


/  ^1 


/  ^1 


8  = 


\^n/ 


I 
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Then  i^ie  note  that  equation  (17)  is  in  the  form: 


(18) 


a\l;  =  G 


where  \|/  and  G  are  vectors  with  M'il  elements,  namely  T^r   .  s.  and 

[    J  >  •'•J 

/g  .  .')■    and  a  is  a  huge  M«N  by  ¥i-d   matrix.   'e  may  write  the 
I  il » ■'-  J 

system  more  compactly  in  tho  form 


(18a) 


Ac|)  =  g 


where    A    is    an   M  by  ¥■  rmtrix  '-nth  elements    that   are   W   by  N 
matrices , 

/'  B      I      0      0    .    .     .    0^.. 


A   = 


I      B      I      0    ...    0 
0      I      B      I    ...    0 


0      0    . 


TBI/ 


0      0    ...    0      I      B 


/ 


where  I  is  the  i^l  by  N  identity  matrix  and  b  is  the  N  by  N  matrix 
of  numbers . 


B  = 


1 

P 
1 

0 
0 


0 

1 

p 


0 
0 

1 


1 

c 


p 
1 


0   \ 
0    \ 

1 

p 


1  with  p 


-ik   +  K^h^) 


Equation  (l8)  always  has  a  unique  solution  regardless  of 
the  value  of  K..  ^   To  show  this  we  consider  the  homogeneous 
equations  a\]/  =  0 .   This  is  equivalent  to  equation  (17)  with 
P.  .  E  0,  and  _T .  .  =  0  on  the  boundary.   Now  this  equation  cer- 
tainly  has  the  trivial  solution.   Suppose  it  had  also  a 


Ik 

non-trivial  solution.   Then  the  value  of  |^ .  .|  for  any  interior 
point  is  not  greater  than  the  average  value  of  |_^ .  .|  for  its 
four  neighboring  points,   '-^e  can  see  this  from  equation  (1?)  and 
noting  -P/U  >  1.   Hence  |JJ.  .  |  is  not  a  positive  maximum  for  any 
interior  point.   This  contradicts  the  assiomption  of  the  exist- 
ence of  a  non-trivial  solution.   Hence  the  homogeneous  equation 
a\l/  =  0  has  only  the  trivial  solution.   This  establishes  the 
unique  existence  of  the  solution  of  (l8). 

To  solve  the  system  (l3)  we  use  the  following  iteration 
method.      First  we  make  a  e'uecs  for  the  values  of  t  •  •  •   In 
practice  this  guess  is  the  value  of  _[,.  .  from  the  previous  time 
step.   Using  superscripts  now  to  denote  the  number  of  the  ite- 

TpTQ 

rant  we  call  this  initial  guess  J.  •  ^  •   Now  we  stort  at  j  =  1 
(the  first  column)  and  solve  the  system: 


(19)  c^S  ^  "^^1  +  ^2  "  Si 


for  the  values 


/iffl    \ 

^i  =   I     :        I      , 
\    *i      ' 

using  the  known  boundary  values  4^  and  the  known  (previously 
guessed)  values  for  <^^.      'e  call  the  solution  (p,  the  first  ite- 
rant for  column  1,  i.e.,  (j)..  =  (J),  .   Then  we  proceed  similarly  to 
solve  column  2  using  c|)-,  +  h^p   +  c|)-  =  f^^p.      '''e  solve  for  4p>  using 
the  freshly  computed  known  values  of  (|),  and  the  old  values  4^. 
Proceeding  thus,  we  solve  for  each  column  and  obtain  the  first 
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Iterants  ^..      Then  we  proceed  in  general  to  solve  for  higher 
iterants  in  a  similar  manner  using: 


(20,       -,     4--  .  Bi'J"  >  i^,,  =  Gj 

ere    c|)'.    ,    and.    ct  . . ,    are   known   either   from  previous    calculations 
J-1  J^-1 


wh 

or  are  knovn  boundary  values.   V''e  solve  equation  (20)  in  the 
following  manner:   If  we  take  the  first  and  third  terms  of  the 
left  side  of  (20)  to  the  right  hand  side,  (20)  is  in  the  form: 

(20a)  Bf^-   --   Z'\^^ 

J      J 

with  S^   =  G .  -  (t^  T  -  ci'^.T  •   i'ow  the  matrix  B  may  be  factored 
J      J     J-1     J+1 

Into  two  Tiatrlces,  i.e.,  B  =  LU,  with  L  a  triangular  matrix  with 
zeroes  above  the  principal  diagonal;  and  U  also  a  triangular  ma- 
trix, but  with  zeroes  below  tlie  principal  diagonal.   Then  set- 

.  .    ttIH' 1    ^n+l    , 
ting  Uq)  .    =    K.  •        we  have: 

(20b)  L^T'l  =  S^"^ 

which  is  simple  to  solve  for  ^^f   ,  since  L  is  a  triangular  ma- 
trix.   Then  we  have: 

(20c)  Uctf  1  =  kX^ 


in+1 
matrix.   A  factorization  of  6  is  given  by 


which  is  again  easy  to  solve  for  4'-    since  U  is  a  triangular 
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(20d) 


^1 

1 

0 

0 

0 

a^ 

1 

0 

0 

0 

a-. 

1 

u  = 


\       0         0 

\ 

\  0        0 


/I 

0 

0 

0 

/ 

1 

^1 

1 

0 

0 

/ 

0 

1 

^2 

1 

0 

L   =' 


0         C 


\0         0 


0      a,_3_      1 


'A   / 


%-2 
0 


^N-1 


1     / 


where    the    factors    a      are    n:lven   recurs  ivelv   bv 

n     '-  -   J 


_  o    1 
^n+1  ~  P  "  a 


and   a-,  -  p 


No  difficulty  in  computinp  these  factors  can  be  encountered 
since  the  sequence  ja  I  is  a  irionotonic,  increasing  sequence  with 
bounds : 


(20e) 


-1+  <  a   <  -2 
^    n 


for  all  H   and  n 


To  accelerate  matters  vje  use  the  extrapolated  Liebmann 
method  and  take  for  the  (n-tl)st  iterant  not  the  solution  of 
equation  (20)  but  a  linear  combination  of  the  nth  iterant  and 
the  solution  of  equation  (20),  i.e.,  we  solve: 
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(21) 


.1-1     ^.1      ^.1  +  1 


in  +  1 
J 


for  <^^,*^   and  then  take  4^"" '''  as 


(22) 


J  J  J 


where  a  is  an  extrapolation  constant  to  be  determined  below. 
1  <  a  <  2.   i^iow  we  have  a  sequence  of  iterants  which  if  they 
converpe,  converpe  to  the  solution  of  equation  (l8). 
If  we  write  equation  (l8a)  in  the  form: 


(23) 


with 


a4  =  C(|)  +  Del)  =  g 


C  = 


/B 

0 

0 

• 

• 

0 

/  I 

j3 

0 

• 

• 

0 

f   ° 

I 

B 

• 

• 

0 

\° 

0 

•    • 
• 

»         • 

I 

B 

0 

\0 

0 

« 

1         • 

0 

T 

B., 

and 


D 


0 

T 

0 

•        • 

0 

\ 

0 

C 

I 

.     . 

0 

0 

0 

0 

»        • 

0 

•    • 

0 

0 

• 

.     .    0 

0 

•^  / 

0 

0 

. 

.     .    0 

0 

0    .' 

and  now  write  e^  =  (4-4^'),  where  4"^  is  the  nth  iterant  of  equa- 
tions (21)  and  (22);  we  may  now  write  the  error  committed  in  the 
(n+l)st  iterant  nn  terms  of  the  error  committed  in  the  nth 
iterant  as: 
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(21|) 


e^^l  =  -C-^De"  =  Ee"^ 


Here,  the  matrix  E  is  equal  to  -C"-'-D   . 


given  by 
(25) 


The  absolute   value   of   the  maxiTium  eigenvalue    of  E   is 
2 


2  2 

a     cos      ('k/M  +   1) 


[p  +   2   cos    (ti/N  +   1)  ] 


2    <   1 


for  all  a  in  1  <  a  <  2.   This  value,  and  hence  the  number  of 
iterations  required  for  a  certain  accuracy,  is  minimized  by 
taking  a  as  the  root  of: 

(26)   a^cos^(7i/M  + 1)  -  (a-l)[p  +  2cos  (ti/N+I)]^  =  0  ,   1  <  a  <  2. 

The    inequality)-    (25)    assures   convergence   of  the    iterants   to 
the   solution  of  equation   (l8). 

Experience    in   actual   computation  shows    that   the   norm   of 
the   error   is   reduced,    on  the   avora,-e,   by   about   2>Q  per  cent   at 
each   iteration  —  requiring,    in  practical   calculations    about  20   " 
iterations    to  produce    6   figure   accuracy. 

VI.      The    Ordinary  Differential  Er-uations    (7b) 


cix^/i-^    =    -(ll^y 


(7b) 


^o'^o'' 


dy^dx   .  (1^)^ 


^o'^o'' 


Since  we  have  available  from  the  computations  (?-,)•  ^  and 
(x  ,y  )",  we  find  (x^,y^_)   "^  in  the  following  manner,  h'e   choose 


o '  '^  o  ■ 
the  four  points  of  the  mesh  surrounding  the  point  (x  ,y 


,k 
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Using  centered  differences  of  (I,  )^  .  we  evaluate  (?,  )   and 

-^1  ,) ,  1  -^1  X 

(.1]_)^  for  each  of  these  four  points.   Then  we  linearly  interpo- 

late  these  four  values  to  find  ( f , )   and  ( Vn )   at  the  point 

—1  X     — 1  y 

(x^,y^)  .   ''6  then  compute: 


dx 


^o"  =  'o  ^  oT^  /at  -  ocAc^)  =  xj  -  2^     Ac  *  o(Ar2) 


(27) 


^r'  =  ^0  ^  ^  At  ^  o(At^,  =  ,1  .11  At  -  ccAx^: 

X 


Then  with  this  approximate  position  (x  ,y  )    ,  we  find 


the  boundary  values  Y^   =  -{y/Zti)}-.    (>(  r  )  where 


-0 


r^  =  /  v'^-^'q)   +  (y-Jq)   ''^nd  (x,y)  are  points  on  the  boundary. 
Hence  we  may  determine  the  boundary  values  of  _^,"   from  the  re- 
lation Ij'l  =  |kH-l_|k+l^   ^,^^  boundary  values  of  f^^^   are   used  to 

Tnk+1 
compute  _Y*    at  the  interior  points. 

V'e  now  difference  and  interpolate,  as  above,  ^'   at  the 

four  points  surrounding  (^o»yo)    >  thus  obtaining  _|-|    and 

_f,    at  the  point  (^  ,y  )    .   ''e  then  compute  finally: 


(28) 


■X  y 

VII.      Computation  of   the   Function  J^  =    -( Y/SirjK^lfcr    ) 

The   strength  of   the    vortex,    Yj    is    taken  as    constant, 
although  provisions   have   been  made   to    allow   tnis   to  vary  in  time 
as   the    computation  progresses.      It    nay   be   calculated   from  the 


20 
initial  data  by  evaluating  the  numerical  analogue  of  the 

0 

Integral  formulg. : 

(29)        .      Y  -    I^  ds  -  k^   ^  ^^^y 

•'b  R 

where  R  is  any  region  containing  the  vortex  and  B  is  its  bound- 
ary.  Another  procedure  is  to  choose  that  value  of  y  such  that 
the  resulting  steering  flow  computed  from  ^^    ~  1  -  ^     -^-S  smooth 
and  has  at  the  position  of  the  vortex  its  contour  lines  tangent 
to  the  direction  of  the  initial  velocity  of  the  vortex. 

The  computation  of  the  K  Bessel  function  as  carried  out 

o 

by  an  automatic  computer  is  effected  by  evaluating  the  classical 
series  representation.   Since  the  range  of  the  argument,  *<r  ,  is 
between  0  and  2,  only  eight  terms  of  the  series  are  needed  to  get 
ten  dlpit  accuracy.   The  slncularity  of  the  function  K   for  zero 
argium.ent  is  circi.imvented  by  trunc^^ting  the  function  when  the 
argument  is  close  to  zero.   In  pi-aotice  this  truncation  is  taken 
at  a  height  on  the  graph  of  '<   which  allows  tangential  velocities 
due  to  the  vortex  no  more  than  a  certain  maximum  on  the  order  of 
100  to  200  m.iles  per  hour.   This  maximum  allowable  velocity  is 
computed  from  the  computational  stability  requirement  given  in 
equation  ( 1L|.)  . 

VIII.   Computation s  on  the  Up. i vac 

The  procedures  outlined  here  for  the  numerical  computation 
were  effected  on  the  Unlvac  machine  at  New  York  University.   A 
coded  program  of  considerable  flexibility  exists  now  to  perform 
automatically  the  computations  and  test  the  feasibility  of  this 
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method.   The  code  will  handle  gilds  of  sizes  up  to  20  by  20 
points.   It  Is  an  especially  nice  feature  of  this  problem  that 
the  variation  of  the  physical  quantities  involved  can  be  esti- 
mated a  priori  and  thus  the  problem  has  been  scaled  accordingly 
for  fixed  point  arithmetic. 

In  pr3ctice,  tl'ne  increments  on  the  order  of  one  hour  with 
mesh  spacing  ropresentinp  distances  on  the  order  of  200  miles 
have  been  used.   The  initial  data  required  is  obtained  by  inter- 
polating from  the  U.  S.  '-/eather  Bureau's  charts  of  the  50C  mil- 
libar constant  pressure  surface.   Detailed  information  concern- 
ing physical  measT.irements  in  the  region  of  a  hurricane  have  been, 
at  the  present  writin.r,  eitlier  non-existent  or  not  obtainable  in 
a  readily  useful  form. 

This  problem  was  coded  for  Univac  rather  than  for  a  faster 
IBM  machine  of  greater  internal  memory  storage  because  one  re- 
tains greater  control  of  computational  flow  in  such  a  pilot  study. 
But  even  with  the  relatively  slow  Univac  and  using  the  coded 
program  at  maxi'num  capacity  a  211.  hour  prediction  is  obtainable 
in  about  L|.  to  5  )-ours .   According  to  our  most  reliable  estimates 
the  IBM  70)j.  can,  in  a  production  run,  produce  a  2L|.  hour  predict- 
ion, using  this  m.ethod ,  in  about  i;  to  '^   minutes.   This  is  possi- 
ble for  two  reasons.   First,  the  IBM  computes  ten  tim.es  faster 
and  secondly  a  factor  of  six  is  saved  by  handling  all  data  in- 
ternally in  t]ie  m.emory  rather  than  trans fei-rin,?:  from  tape  to 
memory  —  a  time-consuming  process  required  in  the  Univac  program. 


22 

IX.   Test  Runs 

The  coded  program  was  tested  ap:ain3t  a  simple  exact  solu- 
tion of  the  peostrophic  conservation  equation.   This  served  se- 
veral purposes;  first,  to  test  the  proper  functioning  of  the 
code  in  all  its  parts,  and,  secondly,  to  verify  the  estimates  of 
the  accuracy  of  computation,  and  thirdly,  to  test  the  proper 
scaling  of  physical  quantities.   The  exact  solution  used  is: 

(30)        I(x,y,T)  =  -Ae*'-y  +  BK^(VTr^)        A,B  constants   . 

This  is  already  in  the  form  ^  =  ?-,  +  ?  •   ^'e  can  compute 
from  (7b)  the  trajectory  of  the  vortex  as: 


(31) 


X  (x)    =  X  (0)  +  KkT  =    kAX 
o   ■     o 


V  (t)  =  y  (0)  =  0   . 
""  o       ''o 


We  computed  and  chose  A  so  that  the  vortex  should  move  with 
velocity,  KA ,    20  m.p.h.  due  east.   Using  time  steps  of  one  hour 
and  mesh  width  of  200  miles  for  a  20  by  20  grid  we  inserted  the 
initial  and  boundary  values  appropriate  to  equation  (30).   The 
computational  results  showed  that  the  vortex  did  indeed  move  20 
m.p.h.  due  east  and  the  steering  flow  reraained  unchanged  in  time, 
as  in  equation  {30),      This  calculation  was  carried  out  for  2[|. 
time  steps  and  the  accuracy  of  the  results  remained  well  within 
the  limits  designed.   Variation  in  qu-antities  was  amply  allowed 
for  in  the  scaling  of  the  physical  quantities. 

Next,  data  estimated  from  a  standard  V?eather  Bureau  chart 
was  used  as  initial  and  boundary  data  for  the  code.   The  results 
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of  this  test  run  for  I4.8  hourly  time  steps  are  being  prepared  and 
partial  results  will  be  analyzed  in  a  later  report  by  G.  K. 
Norikawa. 

At  present,  using  the  expeiaence  gained  in  the  previous 
tests,  we  are  making  a  full-scale  test  of  the  code  on  as  accu- 
rate ^'eather  Bureau  data  as  we  can  obtain  —  specifically,  choos- 
ing the  data  from  an  actual  hurricane,  lone  of  1955 •   A  number 
of  auxiliary  routines  are  beinjj;  coded  to  aid  in  the  editing  and 
preparation  of  both  input  and  output  data. 

X.   Concluding  Remarks 

The  numerical  soli.ition  of  the  geostrophic  conservation 
equation  as  outlined  here  is   the  subject  of  a  forthcoming  doctor- 
al dissertation  by  Mrs.  H.  hontvila.   Ghe  discusses,  in  a  rigo- 
rous m.anner,  the  convergence  of  the  numerical  solution  to  the 
actual  solution  of  the  equations,  for  a  wide  class  of  initial 
and  boundary  value  data,   ^'e  have  previously  m.entioned  the  rigo- 
rous existence  and  uniqueness  proofs  of  C.  Sensenig. 

The  computational  results  of  test  rvms  so  far  strongly 
indicate  the  feasibility  in  both  accuracy  and  speed  of  solving 
the  geostrophic  conservation  equation  by  numerical  methods. 
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